Notizbuch 1: Zeitbeiwertverfahren, Einheitsganglinie, Block- und detaillierte Modelle#
Einführung in die Einheitsganglinie, das Zeitbeiwertverfahren und die konzentrierte und verteilte hydrologische Modellierung.
Einführung#
In diesem Notizbuch betrachten wir verschiedene Möglichkeiten zur Schätzung einer Niederschlags-Abfluss-Ganglinie mit unterschiedlichem Komplexitätsgrad.
(Quelle: Margulis 2020)
Zunächst verwenden wir das Zeitbeiwertverfahren, um anhand empirischer Werte den Spitzenfluss für ein kleines Einzugsgebiet (EZG) im Kanadischen Whistler zu schätzen. Anschließend stellen wir einige komplexere Modelle vor, darunter die Modelle Soil Conservation Service (SCS), Unit Hydrograph (SCS-UH) und Curve Number (SCS-CN), welche ihr auch in der Vorlesung kennengelernt habt. Wir berechnen die maximale Abflusspende mithilfe der CN-Methode, wobei davon ausgegangen wird, dass die EZG-Eigenschaften wie Neigung und Rauheit im gesamten EZG homogen sind. In diesem Fall werden also alle Bereiche des Beckens als eine Einheit behandelt, was oft als Block-Modell bezeichnet wird.
Abschließend verwenden wir eine Open-Source-Geodatenbibliothek, um aus digitalen Höhendaten (DEM) ein einfaches detailliertes Modell des EZGs zu erstellen. Wir berechnen die Fließrichtung und die Fließakkumulation, um ein EZG abzugrenzen und das Gewässernetz zu definieren, und erstellen daraus eine Ganglinie für ein Niederschlagsereignis.
In allen Fällen werden wir die resultierende Abflüsse mit einem praktischen Problem in Zusammenhang bringen, das den Wasserstand in einem hypothetischen Kanal am Beckenauslass betrifft. Hierdurch soll ein Gefühl für die Bandbreite der Ergebnisse (den Spitzenfluss der Niederschlags-Abfluss-Reaktions) unter unsicheren Eingabevariablen entwickelt werden.
# importiere die nötigen Pakete
import os
import pandas as pd
import numpy as np
import math
# fortgeschrittene Statistik Bilbiothek
from scipy import stats
import matplotlib.pyplot as plt
import matplotlib.patches as mp
import matplotlib.colors as colors
from pysheds.grid import Grid
working_directory = os.getcwd()
print(working_directory)
C:\Users\GJ\OneDrive - Universitaet Duisburg-Essen\GitHub\repos\Wasserwirtschaft_und_Klimawandel\Inhalt\Notebooks\Notebook_1
Niederschlagsdaten herunterladen und importieren#
Wir können die Anwendung des Pacific Climate Impacts Consortium (PCIS) nutzen, um Wetterdaten in der Gegend von Whistler abzurufen. Für diese Übung verwenden wir historische Wetterdaten der Station des Meteorological Service of Canada (MSC) in Whistler, BC.
Unter Verwendung des bereitgestellten Filters Obervation Frequency scheint es einige Klimastationen mit stündlichen Niederschlagsdaten zu geben:
Wir betrachten eine Station (ID 1048899: Whistler (2014-2022)) als Beispiel, schauen uns aber auch eine andere an, da die PCIS-Webanwendung nahelegt, dass diese Station stündliche Daten hat, dies aber nicht der Fall ist. Das ist ebenfalls für alle anderen Stationen außer einem (925 – grünes Dreieck, rot umkreist) der Fall. Ihr seid immer für eure eigene Datenvalidierung verantwortlich.
In der PCIS-Datenbank sind nur an einem Standort stündliche Daten um Whistler herrum verfügbar und das nur für einen kurzen Zeitraum im Jahr 2005. Diese Daten sind von der von der Wildfire Management Branch (FLNRO-WMB) des Ministeriums für Wälder, Land und Ressourcen (ID 925 ZZ REMSATWX1).
# importiere Niederschlagsdaten
daily_df = pd.read_csv('../../Notebook_Daten/Notebook_1_Daten/Whistler_348_climate.csv', parse_dates=True, index_col=['Date/Time'])
daily_df.dropna(subset=['Total Precip (mm)'], inplace=True)
daily_df.head()
| Station Name | Year | Month | Day | Data Quality | Max Temp (°C) | Max Temp Flag | Min Temp (°C) | Min Temp Flag | Mean Temp (°C) | ... | Total Snow (cm) | Total Snow Flag | Total Precip (mm) | Total Precip Flag | Snow on Grnd (cm) | Snow on Grnd Flag | Dir of Max Gust (10s deg) | Dir of Max Gust Flag | Spd of Max Gust (km/h) | Spd of Max Gust Flag | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Date/Time | |||||||||||||||||||||
| 1976-11-22 | WHISTLER | 1976 | 11 | 22 | NaN | 5.6 | NaN | NaN | M | NaN | ... | 0.0 | NaN | 1.5 | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 1976-11-23 | WHISTLER | 1976 | 11 | 23 | NaN | 5.0 | NaN | 0.6 | NaN | 2.8 | ... | 0.0 | NaN | 12.7 | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 1976-11-24 | WHISTLER | 1976 | 11 | 24 | NaN | 5.0 | NaN | 2.2 | NaN | 3.6 | ... | 0.0 | NaN | 2.5 | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 1976-11-25 | WHISTLER | 1976 | 11 | 25 | NaN | 4.4 | NaN | 0.6 | NaN | 2.5 | ... | 0.0 | NaN | 0.0 | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 1976-11-26 | WHISTLER | 1976 | 11 | 26 | NaN | 2.2 | NaN | -6.7 | NaN | -2.3 | ... | 0.0 | NaN | 0.0 | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
5 rows × 27 columns
# Merke dass die ascii Datei die Zeichenfolge 'None' für NaN
# wir berücksichtigen das während des Datenimports
hourly_df = pd.read_csv('../../Notebook_Daten/Notebook_1_Daten/925.ascii',
header=1, na_values=[' None'], parse_dates=[' time'])
# Merke das das ascii Format die Spaltenüberschriften mit spaces importiert
# die wir bereinigen müssen
hourly_df.columns = [e.strip() for e in hourly_df.columns]
hourly_df.set_index('time', inplace=True, drop=True)
hourly_df.head()
| wind_speed | temperature | wind_direction | relative_humidity | precipitation | |
|---|---|---|---|---|---|
| time | |||||
| 2005-09-28 10:00:00 | 12.9 | 11.6 | 234.0 | 61.0 | 0.0 |
| 2005-09-28 11:00:00 | 14.8 | 11.2 | 229.0 | 68.0 | 0.0 |
| 2005-09-28 12:00:00 | 10.4 | 11.1 | 224.0 | 76.0 | 0.2 |
| 2005-09-28 13:00:00 | 13.8 | 10.5 | 247.0 | 79.0 | 0.0 |
| 2005-09-28 14:00:00 | 13.9 | 10.0 | 245.0 | 80.0 | 0.2 |
Daten plotten#
Hier verwenden wir die Bokeh-Datenvisualisierungsbibliothek, um die Niederschlagsdaten grafisch darzustellen. Die Bibliothek bietet die Flexibilität, Zeitskalen zu vergrößern und zu verkleinern, was bei der Datenexploration und -überprüfung hilft.
from bokeh.plotting import figure, show
from bokeh.io import output_notebook
from bokeh.models import ColorBar, ColumnDataSource
output_notebook()
hourly_source = ColumnDataSource(hourly_df)
daily_source = ColumnDataSource(daily_df)
p = figure(title=f'Niederschlagsdaten', width=750, height=300, x_axis_type='datetime')
p.vbar(x='Date/Time', width=pd.Timedelta(days=1), top='Total Precip (mm)',
bottom=0, source=daily_source, legend_label='Tägl. Niederschlag',
color='royalblue')
p.vbar(x='time', width=pd.Timedelta(hours=1), top='precipitation',
bottom=0, source=hourly_source, legend_label='Stündlicher Niederschlag',
color='gold')
p.legend.location = 'top_left'
p.xaxis.axis_label = 'Datum'
p.yaxis.axis_label = 'Niederschlag [mm]'
p.toolbar_location='above'
show(p)
Use the zoom tool to compare the independent measurements over 28-29 September 2005. Do the daily volumes add up?
Verwenden Sie das Zoom-Tool, um die gemessenen Werte vom 28. bis 29. September 2005 zu vergleichen. Passen die täglichen Volumina hierzu?
hourly_df['day'] = hourly_df.index.day
hourly_tot = hourly_df.groupby('day')['precipitation'].sum()
print(hourly_tot)
daily_check = daily_df['2005-09-28':'2005-09-29']['Total Precip (mm)'].copy()
print(daily_check)
day
28 17.4
29 22.8
Name: precipitation, dtype: float64
Date/Time
2005-09-28 26.2
2005-09-29 17.1
Name: Total Precip (mm), dtype: float64
hourly_df
| wind_speed | temperature | wind_direction | relative_humidity | precipitation | day | |
|---|---|---|---|---|---|---|
| time | ||||||
| 2005-09-28 10:00:00 | 12.9 | 11.6 | 234.0 | 61.0 | 0.0 | 28 |
| 2005-09-28 11:00:00 | 14.8 | 11.2 | 229.0 | 68.0 | 0.0 | 28 |
| 2005-09-28 12:00:00 | 10.4 | 11.1 | 224.0 | 76.0 | 0.2 | 28 |
| 2005-09-28 13:00:00 | 13.8 | 10.5 | 247.0 | 79.0 | 0.0 | 28 |
| 2005-09-28 14:00:00 | 13.9 | 10.0 | 245.0 | 80.0 | 0.2 | 28 |
| 2005-09-28 15:00:00 | 12.1 | 8.9 | 231.0 | 84.0 | 1.0 | 28 |
| 2005-09-28 16:00:00 | 8.9 | 8.5 | 227.0 | 88.0 | 1.2 | 28 |
| 2005-09-28 17:00:00 | 11.3 | 8.2 | 226.0 | 92.0 | 2.2 | 28 |
| 2005-09-28 18:00:00 | 12.3 | 7.9 | 225.0 | 93.0 | 2.4 | 28 |
| 2005-09-28 19:00:00 | 7.4 | 7.8 | 180.0 | 95.0 | 1.6 | 28 |
| 2005-09-28 20:00:00 | 6.2 | 7.7 | 183.0 | 96.0 | 2.0 | 28 |
| 2005-09-28 21:00:00 | 7.6 | 7.9 | 178.0 | 95.0 | 2.6 | 28 |
| 2005-09-28 22:00:00 | 7.1 | 8.0 | 200.0 | 95.0 | 2.6 | 28 |
| 2005-09-28 23:00:00 | 9.2 | 8.1 | 195.0 | 96.0 | 1.4 | 28 |
| 2005-09-29 00:00:00 | 7.7 | 8.2 | 211.0 | 97.0 | 2.4 | 29 |
| 2005-09-29 01:00:00 | 8.4 | 8.4 | 214.0 | 96.0 | 1.2 | 29 |
| 2005-09-29 02:00:00 | 7.7 | 8.6 | 218.0 | 98.0 | 0.8 | 29 |
| 2005-09-29 03:00:00 | 6.5 | 8.6 | 186.0 | 98.0 | 0.4 | 29 |
| 2005-09-29 04:00:00 | 8.3 | 8.5 | 199.0 | 98.0 | 0.8 | 29 |
| 2005-09-29 05:00:00 | 8.6 | 8.6 | 188.0 | 98.0 | 0.6 | 29 |
| 2005-09-29 06:00:00 | 7.6 | 8.6 | 202.0 | 98.0 | 1.2 | 29 |
| 2005-09-29 07:00:00 | 8.0 | 8.6 | 157.0 | 98.0 | 1.6 | 29 |
| 2005-09-29 08:00:00 | 6.2 | 9.0 | 157.0 | 97.0 | 1.6 | 29 |
| 2005-09-29 09:00:00 | 11.3 | 9.2 | 185.0 | 96.0 | 2.0 | 29 |
| 2005-09-29 10:00:00 | 3.6 | 9.2 | 174.0 | 98.0 | 5.4 | 29 |
| 2005-09-29 11:00:00 | 9.3 | 9.6 | 203.0 | 96.0 | 2.8 | 29 |
| 2005-09-29 12:00:00 | 9.8 | 10.1 | 216.0 | 96.0 | 1.2 | 29 |
| 2005-09-29 13:00:00 | 12.9 | 10.6 | 224.0 | 96.0 | 0.8 | 29 |
Der Niederschlagssumme für beide Tage liegt sehr nahe beieinander, aber die Mengen der einzelnen Tage stimmen nicht überein, obwohl diese Stationen nahe beieinander liegen (die Whistler-Station ist das rote Dreieck direkt nördlich des grünen Dreiecks in der Karte oben).
Wenn wir uns zufällig für diesen bestimmten Tag interessieren, welche anderen Informationen könnten in Bezug auf den Energieeintrag relevant sein? Schaut euch die Datenspalten an:
daily_df.columns
Index(['Station Name', 'Year', 'Month', 'Day', 'Data Quality', 'Max Temp (°C)',
'Max Temp Flag', 'Min Temp (°C)', 'Min Temp Flag', 'Mean Temp (°C)',
'Mean Temp Flag', 'Heat Deg Days (°C)', 'Heat Deg Days Flag',
'Cool Deg Days (°C)', 'Cool Deg Days Flag', 'Total Rain (mm)',
'Total Rain Flag', 'Total Snow (cm)', 'Total Snow Flag',
'Total Precip (mm)', 'Total Precip Flag', 'Snow on Grnd (cm)',
'Snow on Grnd Flag', 'Dir of Max Gust (10s deg)',
'Dir of Max Gust Flag', 'Spd of Max Gust (km/h)',
'Spd of Max Gust Flag'],
dtype='object')
Tagestemperatur und dem Schnee am Boden? Wie wäre es mit Flags zur Datenqualität?
climate_check = daily_df['2005-09-28':'2005-09-29'][['Snow on Grnd (cm)', 'Mean Temp (°C)']].copy()
climate_check
| Snow on Grnd (cm) | Mean Temp (°C) | |
|---|---|---|
| Date/Time | ||
| 2005-09-28 | 0.0 | 7.0 |
| 2005-09-29 | 0.0 | 10.1 |
flag_check = daily_df['2005-09-28':'2005-09-29'][['Snow on Grnd Flag', 'Mean Temp Flag', 'Total Precip Flag']].copy()
flag_check
| Snow on Grnd Flag | Mean Temp Flag | Total Precip Flag | |
|---|---|---|---|
| Date/Time | |||
| 2005-09-28 | NaN | NaN | NaN |
| 2005-09-29 | NaN | NaN | NaN |
Keine besondere Kennzeichnung der Daten und kein Schnee auf dem Boden, daher haben wir möglicherweise etwas mehr Vertrauen in die Tagesdaten, obwohl es schwierig ist, das Timing zwischen den beiden Datensätzen in Einklang zu bringen. Es könnte sich um ein Zeitzonenproblem handeln, da einige Behördern die Zeiten in UTC melden und somit 7–8 Stunden voraus sind. Es könnte auch an der Art und Weise liegen, wie die Zeitstempel bei der Aggregation von Daten gehandhabt werden.
Beschreibung des Projektvorhabens#
Stellt euch vor, ihr arbeitet an einem Projekt zur Planierung und Neupflasterung rund um den Parkplatz 4. Aufgabe des Porjektes ist zudem die Installation einer Entwässerung, um das Wasser vom Parkplatz aufzufangen und es in ein Regenrückhaltebecken (RRB) umzuleiten, anstatt es in den FitzSimmons Bach abzuleiten. Hier erfolgt ebenfalls eine stoffliche Vorbehandlung, um bei Einleitung nicht das dortige Ökosystem zu schädigen.
Projektstart ist im Sommer und das Projekt soll im Herbst abgeschlossen sein. Nehmt für diese Aufgabe an, dass es sich um ein Einzugsgebiet von ungefähr 1 \(km^2\) handelt, welches durch einen temporäres Gerinne (zu planen) in ein RRB mündet, bevor es in den FitzSimmons Bach weitergeleitet wird.
Frage: Habt ihr angesichts der Tatsache, dass das Einzugsgebiet nur \(1 km^2\) ist, eine Idee, welche Niederschlagsdauer für die Schätzung Scheitelabflusses erforderlich ist? d.h. 1 Minute, 10 Minuten, 1h, 6h, 24h, 48h?
Es gibt eine Stelle, an der ein natürliches Gerinne vorhanden ist, welches für das Projekt genutzt werden könnte. Hierdurch könnten Ausgrabungen vermieden und der Projektzeitplan beschleunigt werden. Der Standort müsste also vermessen werden, um genaue Abmessungen zu erhalten.
Angenommen, ihr beuftragt die Vermessung und erhaltet von der Vermesser:in folgende Ergebnisse der natürlichen Geometrie.
*Kanalbreite (b): \(3m\)
*Kanaltiefe (h): \(0.5m\)
*Hydraulisches Gefälle (S): 0,005 (0,5 %)
*Rauheit (n): 0,017 (rauer Asphalt)
Berechnen der Abflusskapazität#
Da die Gerinnegeometrie bekannt ist ist es möglich, die Abflusskapazität unseres Grinnes zu bestimmen. Wir berechnen zunächst die maximale Kapazität und bewerten dann, wie sich unser Spitzenabfluss hierzu verhält. Zur Abschätzung der Fließgeschwindigkeit eines Gewässers können wir die Fließformel nach Gauckler-Manning-Strickler anwenden:
b, h = 2, 0.5 # height, width in metres
# Querschnittsfläche
A = b * h
perim = 2 * b + 2 * h
# hydraulischer Radius
R = A / perim
# Gefälle
S = 0.005
# Strickler Beiwert
n_manning = 0.017
# average velocity
V_avg_Gerinne = R**(2.0/3) * S**(1.0/2) / n_manning
print(f'Die durchschnittliche Fließgeschwindigkeit bei vollem Durchfluss beträgt {V_avg_Gerinne:.1f}m/s')
Q_Gerinne_Kapazität = V_avg_Gerinne * A
print(f'Die volle Abflusskapazität des Gerinnes beträgt {Q_Gerinne_Kapazität:.1f} m^3/s')
Die durchschnittliche Fließgeschwindigkeit bei vollem Durchfluss beträgt 1.4m/s
Die volle Abflusskapazität des Gerinnes beträgt 1.4 m^3/s
Berechnung der Ganglinie#
Langzeitaufzeichnungen mit hoher zeitlicher Auflösung sind nur in seltenen Fällen zu finden, weshalb wir mit den verfügbaren Informationen unser Bestes geben. Im Folgenden sehen wir uns einige Möglichkeiten zur Erstellung einer Ganglinie mit unterschiedlichem Detaillierungsgrad an. Grundsätzlich wollen eine Ganglinie erstellen oder zumindest ihren Scheitelabfluss genau vorhersagen, um Wasserbauwerke und andere Wassermanagementsysteme planen zu können. Wir beginnen mit einer sehr einfachen Schätzung, die nur minimale Eingangsinformationen erfordert, und gehen dann zu komplexeren und Eingangsdaten intensiveren Methoden über.
In der Hydrologie interessieren v.a. Unter- und Überschreitungswahrscheinlichkeiten. Mit anderen Worten: Wie groß ist die Wahrscheinlichkeit, dass der Durchfluss an einem bestimmten Ort in einem bestimmten Jahr die Kapazität des Gerinnes (oder der gebauten Struktur) überschreitet? Für diese Art von Problemen gibt es keine richtige Antwort, sie sind offen und subjektiv – was bedeutet, dass jede Antwort ein gewisses Maß an (technischem) Urteilsvermögen erfordert. Das Thema Risiko wird in der Vorlesung besprochen. Im Moment konzentrieren wir uns nur auf verschiedene Verfahren zur Übertragung einer effektiven Niederschlagsganglinie in eine Direktabflussganglinie.
Konvertiert das Volumen (mm pro Stunde oder Tag) in die Einheit des Volumenstromes#
Der Abfluss wird normalerweise in \(\frac{m^3}{s}\) gemessen und der Niederschlag wird häufig in mm pro Stunde oder Tag angegeben. Konvertiert den \(\frac{mm}{Tag}\)-Niederschlag in \(\frac{m^3}{s}\).
Zeitbeiwertverfahren#
Das Zeitbeiwertverfahren ist das am häufigsten verwendete Fließzeitverfahren zur Ermittlung des Scheitelabflusses von einer gleichmäßig überregneten Fläche mit einer konstanten Regenspende und einem konstanten Abflussbeiwert. Im englischsprachigen Raum wird das Verfahren als “rational formula” oder “rational method” (Verhältnismethode) bezeichnet.
Die Bezeichnung “Zeitbeiwertverfahren” entstammt aus früheren Zeiten, als die Regenspende I aus dem Produkt einer Bezugsregenspende und eines Zeitbeiwertes C für eine bestimmte Regendauer und Regenhäufigkeit berechnet wurde.
Info zu dem Vorgehen in Deutschland: Inzwischen werden meistens die regionalisierten Starkniederschlagsauswertungen des Deutschen Wetterdienstes (Kostra-Atlas) oder andere Niederschlagsauswertungen aus örtlich gewonnenen Daten verwendet.
Bei den Fließzeitverfahren wird davon ausgegangen, dass die maßgebende Regendauer der Fließzeit im betrachteten Einzugsgebiet entspricht. Diese Zeit wird auch als Konzentrationszeit bezeichnet und kann auf dieser Seite berechnet werden:
mit:
Q: Scheitelabfluss [\(m^3/s\)]
k: 0.278 [-] Notiz: \(\quad 1\frac{mm}{hr} \cdot \frac{1\text{hr}}{3600s} \cdot \frac{1m}{1000mm} \cdot 1 \text{km}^2 \cdot \frac{1\times 10^6 m^2}{1\text{km}^2} = \frac{1}{3.6} = 0.278 \frac{m^3}{s}\)
C: Abflussbeiwert [-]
I: Regenspende [mm/hr]
a: Einzugsgebietsfläche [\(km^2\)]
Wir haben das Einzugsgebiet bereits ermittelt und benötigen somit nur zwei weitere Informationen, um den Spitzenabfluss für unser Einzugsgebiet schätzen zu können. Das Verhältnis zwischen Abfluss und Niederschlag (der Abflusskoeffizient) “C” ist in einer Tabelle mit empirischen Werten im obigen USACE-Link zu finden, und die Niederschlagsintensität kann anhand der von Environment Canada entwickelten Intensitäts-Dauer-Häufigkeits-Kurve geschätzt werden. Weitere Informationen zur Verwendung von IDF-Kurven findet ihr hier. Wir können IDF-Kurven für bestimmte Orte mit der Webanwendung auf climatedata.ca finden. Die IDF-Kurve für Whistler findert ihr hier:
Wir kennen die richtige Dauerstufe (X-Achse) für unser Einzugsgebiet noch nicht wirklich, aber wir können mehrere Dauerstufen auswählen und eine Sensitivitäsanalyse durchführen, um zu sehen, wie empfindlich unser Modell auf die Dauerstufen reagiert. Die fünf diagonalen Linien stellen unterschiedliche Wiederkehrintervalle (auch Jährlichkeiten) dar (2, 5, 10, 20, 50, 100 Jahre). Das Wiederkehrintervall ist die Umkehrung der Überschreitungswahrscheinlichkeit und stellt wiederum die Eintrittswahrscheinlichkeit in einem bestimmten Jahr dar. ACHTUNG: Das bedeutet nicht, dass ein Ereignis irgendeiner Größenordnung einmal in diesem Wiederkehrintervall auftritt!
Im Folgenden werden wir die Werte, bei denen wir uns nicht ganz sicher sind, als plausible Bereiche ausdrücken und ihre Auswirkung auf das Ergebnis testen.
def Zeitbeitwert_Spitzenabfluss(C, I, a):
"""Berechne den Spitzenabfluss über das Zeitbeiwertverfahren.
Args:
C (float): Abflussbeiwert [-]
I (float): Regenspende [mm/hr]
a (float): Einzugsgebietsfläche [km^2]
"""
return C * I * a / 3.6
# für jedes Wiederkehrintervall werden wir die minimale und maximale Intensität ablesen
# und verwenden, um die Bandbreite der Ergebnisse zu visualisieren
# Dauer (Minuten): Volumen (mm)
IDF_dict = {
5: (22, 45), # die 2- und 100-Jahres-Intensitäten betragen 22 mm bzw. 45 mm / h
15: (14, 28),
60: (7, 15),
1440: (2, 4) # 1440 Minuten sind 24 Stunden
}
# hier definieren wir ein Array mit drei Werten für den Abflusskoeffizienten.
# um ein Gefühl für den Bereich der möglichen Bedingungen zu bekommen
# Der array erhält Minimum, Maximum und erwarteter Wert
C_values = [0.5, 0.7, 0.9]
# Berechnen Sie den Bereich der Abflussschätzungen für jedes C und jede Wiederkehrperiode.
# und erstellen Sie ein Diagramm für jeden Abflusskoeffizienten
figs = []
rational_results = {}
colors=['green', 'orange', 'red']
for k, (i_min, i_max) in IDF_dict.items():
i = 0
p = figure(title=f'Dauerstufe={k}min', width=600, height=400)
for c in C_values:
Q_min = Zeitbeitwert_Spitzenabfluss(c, i_min, 1)
Q_max = Zeitbeitwert_Spitzenabfluss(c, i_max, 1)
print(f'{k} Dauer Minuten, C={c} Abflussbereich = {Q_min:.1f} to {Q_max:.1f} m³/s')
x = [2, 100]
y = [Q_min, Q_max]
p.line(x, y, legend_label=f'C={c}', color=colors[i])
p.yaxis.axis_label = 'Abfluss [m³/s]'
p.xaxis.axis_label = 'Wiederkehrintervall [1/a]'
p.legend.location = 'top_left'
i += 1
figs.append(p)
5 Dauer Minuten, C=0.5 Abflussbereich = 3.1 to 6.2 m³/s
5 Dauer Minuten, C=0.7 Abflussbereich = 4.3 to 8.7 m³/s
5 Dauer Minuten, C=0.9 Abflussbereich = 5.5 to 11.2 m³/s
15 Dauer Minuten, C=0.5 Abflussbereich = 1.9 to 3.9 m³/s
15 Dauer Minuten, C=0.7 Abflussbereich = 2.7 to 5.4 m³/s
15 Dauer Minuten, C=0.9 Abflussbereich = 3.5 to 7.0 m³/s
60 Dauer Minuten, C=0.5 Abflussbereich = 1.0 to 2.1 m³/s
60 Dauer Minuten, C=0.7 Abflussbereich = 1.4 to 2.9 m³/s
60 Dauer Minuten, C=0.9 Abflussbereich = 1.8 to 3.8 m³/s
1440 Dauer Minuten, C=0.5 Abflussbereich = 0.3 to 0.6 m³/s
1440 Dauer Minuten, C=0.7 Abflussbereich = 0.4 to 0.8 m³/s
1440 Dauer Minuten, C=0.9 Abflussbereich = 0.5 to 1.0 m³/s
from bokeh.layouts import gridplot
layout = gridplot(figs, ncols=2, width=350, height=300)
show(layout)
Frage: Ist es nach den obigen Diagrammen wichtiger, die richtige Dauerstufe zu nutzen, ein geeignetes Wiederkehrintervall zu wählen oder einen geeignetetn Abflussbeiwert zu wählen? Mit anderen Worten: Wie empfindlich reagiert das Modell auf die verschiedenen Eingabeparameter?
Quelle für Erläuterungen zu SCS: https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&ved=2ahUKEwjV58zCnYH_AhURX_EDHd5YBjcQFnoECAsQAQ&url=https%3A%2F%2Fabstracts.boku.ac.at%2Fdownload.php%3Fdataset_id%3D18982%26property_id%3D107&usg=AOvVaw06VzVfxdy25ULOW27durCk
SCS Einheitsganglinienverfahren#
Das folgende Einheitsgangliniendiagramm beschreibt eine idealisierte Einzugsgebietsreaktion auf ein Niederschlagsereignis und wird in physikalisch basierten Modellen verwendet, um den Scheitelabfluss abzuschätzen. Die Form der Ganglinie hat viele Einflüsse und wird durch komplexe Wechselwirkungen zwischen der Atmosphäre, der Oberfläche (d. h. Vegetation, Rauheit, Neigung) und dem Untergrund (d. h. vorhergehende Bodenfeuchtigkeit, Durchlässigkeit, Porosität) beeinflusst. Im Allgemeinen werden diese Werte sorgfältig durch gleichzeitige Abfluss- und Niederschlagsmessung kalibriert, die ein möglichst breites Spektrum an Bedingungen abdecken.
Der Informationsbedarf ist bei der SCS-UH-Methode im Vergleich zum Zeitbeiwertverfahren hoch und zudem stark ortsabhängig und zeitlich sehr variabel. Detaillierte Informationen sind oft nur bei Anwendungen wie der städtischen Hydrologie verfügbar, wo detaillierte Vermessungsinformationen, Oberflächenbeschaffenheit, Neigungen usw. erfasst werden. In einem solchen Fall kann das Soil Conservation Service (SCS) Unit Hydrograph Model nützlich sein.
\(U_P\): Spitzenabfluss
\(C\): Peak Rate Faktor
\(A\): Einzugsgebietsfläche
\(T_p\): Zeit vom Beginn des Niederschlages bis zur Abflussspitze
Schematische Darstellung des dimensionslosen SCS-UH. (Quelle:NRCS)
Der maßgebende Parameter, der bei der Anwendung des SCS-UH Verfahrens geschätzt werden muss, ist die Konzentrationszeit tc. Die anderen Parameter ergeben sich aus Eigenschaften des Einzugsgebietes und des Niederschlages, die im Anwendungsfall bekannt sind, oder, wie zum Beispiel der Effektivniederschlag, in vorhergehenden Modellierungsschritten errechnet werden.
Die Konzentrationszeit (engl.: time of concentration) \(t_c\) entspricht der Fließzeit des abflusswirksamen Niederschlags zwischen der Beobachtungsstelle und dem am weitesten entfernten Punkt des Einzugsgebietes.
Im Folgenden wird die Geschwindigkeitsbasierte Methode (Velocity Method) erläutert, welche vom NRCS in Kapitel 15 des National Engineering Handbook 630 empfohlen wird, um die Konzentrationszeit zu schätzen
Es werden dabei drei elementare Fließprozesse unterschieden die in einer zeitlichen Abfolge hintereinander stattfinden: Zunächst erfolgt der Abfluss als verteilter Flächenabfluss („sheet flow“), welcher in den konzentrierten Flächenabfluss übergeht („shallow concentrated flow“) und zuletzt als Gerinneabfluss („open channel flow“) dem Gebietsauslass zustrebt. In der Geschwindigkeitsbasierten Methode wird zunächst für jeden dieser Abflussprozesse die Fließzeit ermittelt. Die Summe der einzelnen Fließzeiten ergibt sodann die Konzentrationszeit \(t_c\):
Es ist zu beachten, dass die Konzentrationszeit wesentlich von den Eigenschaften des Niederschlagsereignisses abhängt. Der NRCS schlägt vor, der Berechnung der Konzentrationszeit das 2-jährliche, 24-stündige Niederschlagsereignis zugrunde zu legen, da davon ausgegangen wird, dass dieses repräsentativ für eine weite Bandbreite von Ereignissen ist.
Zur Ermittlung des verteilten Flächenabflusses wird eine von Welle und Woodward (1986) entwickelte vereinfachte Version der Fließformel nach Manning vorgeschlagen:
mit:
\(N\): empirischer Rauhigkeitskoeffizient [-]
\(L\): Fließlänge des verteilten Flächenabflusses [Fuß]
\(P\): Niederschlagsmenge eines 2-jährlichen, 24-stündigen Ereignisses [Zoll]
\(S\): Oberflächengefälle
Der verteilte Oberflächenabfluss geht in den konzentrierten Flächenabfluss über, der entlang von Rinnen und ähnlichen Vertiefungen an der Oberfläche abfließt, ohne dass jedoch ein durchgehendes Gerinne mit einem abgrenzbaren Querschnitt definiert werden kann. Für die Fließlänge des konzentrierten Flächenabflusses kann folgende Formel verwendet werden.
mit:
\(L\): Fließlänge des verteilten Flächenabflusses [Fuß]
\(S\): Oberflächengefälle
\(V\): Durchschnittliche Fließgeschwindigkeit [ft/s]
Aus Abbildung 3-1 der SCS Curve Number Method documentation können wir die durchschnittliche Fließgeschwindigkeit für einen flachen, konzentrierten Fluss über unbefestigte Oberflächen (Kies) auf 0,37 m/s (1,2 ft/s) schätzen. Bei einer Gerinnelänge von 250 m entspricht dies 0,2 Stunden. Wir können dies nutzen, um die Konzentrationszeit für unser Becken zu schätzen.
Anmerkung: Seid vorsichtig mit den Einheiten! Viele ältere technische Unterlagen sind in imperialen Einheiten (ft-lb) geschrieben. Wir müssen sicher sein, dass empirische Koeffizienten keine Fehler in unsere Berechnungen bringen.
def calculate_t_sheet(n, L, P, s):
"""
n: empirischer Rauhigkeitskoeffizient [-]
s: Gefälle (ft/ft)
L: Fließlänge (ft)
P: Niederschlagsmenge eines 2-jährlichen, 24-stündigen Ereignisses [Zoll]
"""
return 0.007 * (n * L)**0.8 / ((P**0.5) * s**0.4 )
t_t = 0.2 # 0.2 Stunden Fließzeit
n = 0.0011
s = 0.005
L = 330 # 100m entsprechen ~ 100 ft.
# Aus der IDF-Kurve ergibt sich ein für ein 2-jähriges Ereignis ein 24h-Volumen von 2mm/h * 24h = 48mm = ~1,9 Zoll
P_in = 1.9
t_sheet = calculate_t_sheet(n, L, P_in, s)
print(f'Die Konzentrationszeit ist ~ {t_sheet+t_t:.2f} Stunden: {t_sheet:.2f}h (sheet) + {t_t:.2f}h (shallow).')
Die Konzentrationszeit ist ~ 0.22 Stunden: 0.02h (sheet) + 0.20h (shallow).
Above, we estimated the time of concentration for this very small basin which should help reduce the range of peak flow estimates calculated above, where we assumed a range of 5 minutes to 24 hours duration to get our values for intensity. We might be justified in narrowing our range between 10 and 30 minutes or so if we are confident in the time of concentration estimate.
Oben haben wir die Konzentrationszeit für unser sehr kleines Einzugsgebiet ermittelt. Basierend auf den Werten, könnten wir unser Dictionary der verschiedenen Dauern (5 Minuten bis 24 h) bei dem oben berechneten Spitzenflusss verringern. Es könnte also sinnvoll sein, diesen Bereich auf etwa 10 bis 30 Minuten einzuschränken, wenn wir uns auf die Ermittlung der Konzentrationszeit verlassen können.
SCS Curve Number (CN)-Verfahren#
Beim SCS Curve Number Verfahren wird die erforderliche Datengrundlage ebenfalls gering gehalten. Es verfügt aber im Vergleich zum Zeitbeiwertverfahren über mehr Parameter durch die die hydrologischen Bodeneigenschaften (Bodenart, Bodenfeuchte, Bodennutzung, Jahreszeit) zur Schätzung von Oberflächenverlusten parametrisiert werden können (siehe Kapitel 2 im verlinkten PDF):
\(P\): Niederschlagsmenge [in/[Zeit]]
\(S\): potenzieller Speicher nach Abflussbeginn [in]
\(I_a\): Anfangsverluste vor Abflussbeginn [in] (\(I_a \approx 0.2S\))
CN: Kurvennummer CN („Curve Number“)
Aus entsprechenden Tabellenwerken geht hervor, dass CN-Werte für Schotterstraßen für die hydrologischen Bodengruppen A, B, C und D zwischen 76 und 91 liegen.
Wir müssen auch die Konzentrationszeit abschätzen, können hierfür aber unsere zuvor berechnete Konzentrationszeit verwenden.
def CN_methode_Q(p_mm, cn):
"""Berechne Spitzenabfluss Q (m³/s) basierend auf der CN.
CN is die dem Bodentypen entsprechende curve number
S ist ein Faktor für den Gebietsspeicher nach Abflussbeginn
P ist eine Niederschlagsmenge
"""
S = (1000.0 / cn) - 10
p_in = p_mm / 25.4 # Umrechnung in inches
V_in = (p_in - 0.2 * S)**2 / (p_in + 0.8 * S) # in/h
# 1 inch = 25.4 mm
V_mm = V_in * 25.4 # mm/h
# Umrechnung von mm to m³/s für 1 km²
V_cms = (V_mm / 3600) * (1E6 / 1000)
return V_cms
CNs = [76, 85, 89, 91]
# Die Niederschlagswerte werden als Bereich geschätzt.
# aus den obigen IDF-Kurven für 15 min und 60 min
# Dauer für eine Intensität mit einer Wiederkehrperiode von 2 bis 100 Jahren
p_schätzungen_mm = [12, 15, 28] # mm/h
p_bereich = list(range(15, 50, 5))
CN_Q = {}
for p in p_bereich:
CN_Q[p] = []
for cn in CNs:
V_cms = CN_methode_Q(p, cn)
CN_Q[p].append(V_cms)
Darstellung der Ergebnisse des CN-Verfahrens#
cn_df = pd.DataFrame(CN_Q)
fig, ax = plt.subplots(figsize=(6, 4))
for p in p_bereich:
ax.plot(CNs, cn_df[p], label=f'{p}mm')
ax.set_xlabel('CN')
ax.set_ylabel('Spitzenabfluss [m³/s]')
ax.legend()
<matplotlib.legend.Legend at 0x10e077fb100>
Flächendetailliertes Modell aus räumlichen Daten#
Wie in der Vorlesung besprochen, braucht der Niederschlag Zeit, um zum Auslass eines EZGs zu fließen. Auf dem Fließweg gibt es komplexe Wechselwirkungen mit der Atmosphäre, der Vegetation, der Bodenoberfläche und dem Untergrund, die es sehr schwierig machen, die Reaktion des Einzugsgebiets auf Niederschläge genau vorherzusagen. Jetzt verwenden wir DEM-Daten, um ein Becken abzugrenzen und ein einfaches flächendetailliertes Modell zu erstellen, um aus stündlichen und täglichen Niederschlagsdaten eine Einheitsganglinie zu erstellen. Wir versuchen im Anschluss die von unserem Modell vorhergesagte Abflussmenge mit einer Zeitreihe gemessener Werte aus dem Fitzsimmons Creek zu vergleichen.
Hinweis: Von hier an verlassen wir das Beispiel des 1 \(\text{km}^2\)-Einzugsgebietes und betrachten ein größeres EZG der gleichen Region, das jedoch in den Fitzsimmons Creek entwässert.
Importieren und visualisieren des DEM#
Hier verwenden wir die Bibliothek Pysheds, um das DEM zu verarbeiten.
data_path = '../../Notebook_Daten/Notebook_1_Daten/'
dem_file = 'Whistler_DEM.tif'
dem_path = data_path + dem_file
grid = Grid.from_raster(dem_path)
dem = grid.read_raster(dem_path)
fig, ax = plt.subplots(figsize=(8,6))
fig.patch.set_alpha(0)
plt.imshow(dem, extent=grid.extent, cmap='terrain', zorder=2)
plt.colorbar(label='Höhe (m)')
# Zeige den Ort Whistler als Punkt auf dem Einzugsgebiet
plt.scatter(x=[-122.9535], y=[50.1162], zorder=3, c='red', s=40, label='Whistler Tal')
plt.grid(zorder=0)
plt.title('Digitales Höhenmodell', size=14)
plt.xlabel('längengrad')
plt.ylabel('Breitengrad')
plt.legend()
plt.tight_layout()
Identifizierung und Füllung von Senken im DEM#
Die bevorzugte Eingabe für den Fließrichtungs-Prozess ist ein digitales Höhenmodell (DEM) ohne Senken. Wenn Senken vorkommen, kann dies zu einem fehlerhaften Fließrichtungs-Raster führen. In einigen Fällen können die Daten identifizierte Senken enthalten. Deswegen ist es i.d.R. wichtig, die Morphologie der Fläche gut genug zu kennen, um zu wissen, welche Features wirklich Senken in der Erdoberfläche sind, und bei welchen Features es sich nur um Fehler in den Daten handelt
# Aufbereitung DEM
# ----------------------
# Senken auffüllen
pit_filled_dem = grid.fill_pits(dem)
# Mulden auffüllen
flooded_dem = grid.fill_depressions(pit_filled_dem)
# Auflösen flacher Bereiche für eine kontinuierliche Abwärtsströmung
inflated_dem = grid.resolve_flats(flooded_dem)
Fließrichtungen aus DEM ableiten#
Jeder Zelle eines DEM-Rasters wird eine Fließrichtung zugewiesen, die durch Bestimmung der Neigung an diesem Pixel/dieser Zelle berechnet wird. Die Steigung wird durch Höhenvergleich mit den 8 umgebenden/verbundenen Pixeln ermittelt. Den Richtungen wird einer von acht Werten zugewiesen, die 8-Bit-Ganzzahlen (\(2^1, 2^2, 2^3, \dots 2^7\)) entsprechen.
32 |
64 |
128 |
16 |
- |
1 |
8 |
4 |
2 |
import matplotlib.cm as cmx
import matplotlib.colors as colors
# D8-Flussrichtungen aus DEM bestimmen (D8 = Fließmodellalgorithmus)
# ----------------------
# Richtungsspezifisches Mapping festlegen
dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
# Berechne Fließrichtungen
# -------------------------------------
fdir = grid.flowdir(inflated_dem, dirmap=dirmap)
fig = plt.figure(figsize=(8,6))
fig.patch.set_alpha(0)
# Wir müssen die Farben normalisieren umd Kontrast herzustellen
cNorm = colors.PowerNorm(gamma=0.5)
plt.imshow(fdir, extent=grid.extent, norm=cNorm,
cmap='terrain', zorder=2)
boundaries = ([0] + sorted(list(dirmap)))
plt.colorbar(boundaries= boundaries,
values=sorted(dirmap))
plt.xlabel('Längengrad')
plt.ylabel('Breitengrad')
plt.title('Fließrichtungs Raster', size=14)
plt.grid(zorder=-1)
plt.tight_layout()
Abflussakkumullation ableiten#
Da Jeder Zelle ist eine Richtung zugeordnet ist, ist das Ergebnis der Abflussakkumulation ein Raster zu jeder Zelle, bestimmt durch das Akkumulieren der Gewichtung für alle Zellen, die in einzelne tiefer gelegene Zellen fließen.
# Berechne Abflussakkumullation
# --------------------------
acc = grid.accumulation(fdir, dirmap=dirmap)
fig, ax = plt.subplots(figsize=(8,6))
fig.patch.set_alpha(0)
plt.grid('on', zorder=0)
im = ax.imshow(acc, extent=grid.extent, zorder=2,
cmap='cubehelix',
norm=colors.LogNorm(1, acc.max()),
interpolation='bilinear')
plt.colorbar(im, ax=ax, label='höher gelegene Zellen')
plt.title('Abflussakkumulation', size=14)
plt.xlabel('Längengrad')
plt.ylabel('Breitengrad')
plt.tight_layout()
Abgrenzung eines EZGs über Wasserscheiden#
Nachfolgend skizzieren wir ein Einzugsgebiet, das dem Pegel am Fitzsimmons Creek entspricht. Water Survey Kanada hat kürzlich Polygone von Einzugsgebieten für fast 7000 Pegel in ganz Kanada veröffentlicht. Darüber hinaus werden Pegel auch als shape files bereitgestellt. Eine Einschränkung ist leider, dass wir unseren Punkt als Datei (.shp oder .geojson) angeben müssen und nicht nur Koordinaten bereitstellen können. Der Referenzpegel soll hier der Auslauf des Einzugsgebietes sein. Wenn wir hochaufgelöste Daten und unvollständige Koordinaten eines Abflusspunkts haben, erhalten wir nicht das richtige Pixel, das dem Auslass entspricht.
# Einzugsgebiet abgrenzen
# ---------------------
# Abflusspunkt festlegen (WSC Station am Fitzsimmons Creek in Whistler: 08MG026)
x, y = -122.9488, 50.12025
# Hinweis: Wenn mithilfe des Schwellenwertes ein Einzugsgebiet definiert wird,
# stellen die Abflusspunkte des EZGs die Knoten eines Wasserlaufnetzes dar,
# das aus der Abflussakkumulation abgeleitet wird. Deshalb müssen sowohl ein
# Abflussakkumulations-Raster als auch die Mindestzahl von Zellen, aus denen ein
# Wasserlauf besteht (der Schwellenwert), angegeben werden.
# einen Schwellenwert für die Abflussakkumulation in Zellen angeben
# Für dieses Raster entsprechen 500 Zellen etwa 0,5 km².
accumulation_threshold = 500
# Abflusspunkt zu hoch akkumulierter Zelle fangen (snappen)
x_snap, y_snap = grid.snap_to_mask(acc > 1000, (x, y))
# Einzugsgebiet abgrenzen
catch = grid.catchment(x=x_snap, y=y_snap, fdir=fdir, dirmap=dirmap,
xytype='coordinate')
# Zurechtschneiden und Plotten des Einzugsgebiets
# ---------------------------
# Bounding Box auf das Einzugsgebiet kürzen (clippen)
grid.clip_to(catch)
clipped_catch = grid.view(catch)
# Einzugsgebiet plotten
fig, ax = plt.subplots(figsize=(8,6))
fig.patch.set_alpha(0)
plt.grid('on', zorder=0)
im = ax.imshow(np.where(clipped_catch, clipped_catch, np.nan), extent=grid.extent, zorder=1, cmap='Greys_r')
plt.xlabel('Längengrad')
plt.ylabel('Breitengrad')
plt.title('Abgegrenztes Einzugsgebiet', size=14)
Text(0.5, 1.0, 'Abgegrenztes Einzugsgebiet')
# Flussnetzwerk extrahieren
# ---------------------
branches = grid.extract_river_network(fdir, acc > 500, dirmap=dirmap)
fig, ax = plt.subplots(figsize=(8.5,6.5))
plt.xlim(grid.bbox[0], grid.bbox[2])
plt.ylim(grid.bbox[1], grid.bbox[3])
ax.set_aspect('equal')
for branch in branches['features']:
line = np.asarray(branch['geometry']['coordinates'])
plt.plot(line[:, 0], line[:, 1])
_ = plt.title('D8 Abflussnetzwerk', size=14)
Berechnen der Entfernung zu den flussauffärtigen Zellen#
# Berechnung der Entfernung von jeder Zelle zum Auslass
# -------------------------------------------
dist = grid.distance_to_outlet(x=x_snap, y=y_snap, fdir=fdir, dirmap=dirmap,
xytype='coordinate')
fig, ax = plt.subplots(figsize=(8,6))
fig.patch.set_alpha(0)
plt.grid('on', zorder=0)
im = ax.imshow(dist, extent=grid.extent, zorder=2,
cmap='cubehelix_r')
plt.colorbar(im, ax=ax, label='Distanz zum Auslauf (Zellen)')
plt.xlabel('Längengrad')
plt.ylabel('Breitengrad')
plt.title('Fließlänge', size=14)
Text(0.5, 1.0, 'Fließlänge')
Berechnen der gewichteten Fließzeit#
Basierend auf der Annahme, dass sich Wasser mit derselben Geschwindigkeit fortbewegt (Oberflächenabfluss ist langsamer), bis es sich in einem Gerinne ansammelt und seine Geschwindigkeit startk ansteigt, ordneen wir jeder Zellte eine Fließzeit zu.
# Berechnung der Abflussakkumulation
acc = grid.accumulation(fdir)
# Wir nehmen an, dass Wasser in unseren Gerinnezellen (>= 500 Akkumulation)
# 10 mal schneller fließt als in Hangzellen (< 500 Ansammlungen)
# d.h. wenn die durchschnittliche Fließgeschwindigkeit 1m/s beträgt, beträgt die Hanggeschwindigkeit 0,1m/s = 10 cm/s
weights = acc.copy()
weights[acc >= accumulation_threshold] = 0.1
weights[(acc > 0) & (acc < accumulation_threshold)] = 1
# Berechne gewichtete Entfernung zum Auslass
weighted_dist = grid.distance_to_outlet(x=x_snap, y=y_snap, fdir=fdir, weights=weights, xytype='coordinate')
fig, ax = plt.subplots(figsize=(8,6))
fig.patch.set_alpha(0)
# plt.grid('on', zorder=0)
im = ax.imshow(weighted_dist, zorder=2,
cmap='cubehelix_r')
plt.colorbar(im, ax=ax, label='Distanz zum Auslass (Zellen)')
plt.xlabel('Längengrad')
plt.ylabel('Breitengrad')
plt.title('Gewichtete Distanz zu Auslass', size=14)
Text(0.5, 1.0, 'Gewichtete Distanz zu Auslass')
Entwickeln Sie ein einfaches Niederschlags-Abfluss-Modells#
Da wir nun die gewichteten Fließentfernungen für jede Zelle im abgegrenzten Einzugsgebiet haben, können wir den Niederschlag auf jede „Zelle“ anwenden, um eine Abflussganglinie zu rekonstruieren.
Zuerst benötigen wir die Zellenabmessungen (die Größe jedes Rasterpixels auf dem Boden). Aus den USGS Hydrosheds-Informationen wissen wir, dass die Auflösung 15 (Grad) Sekunden beträgt. Da die Erde keine perfekte Kugel ist, werden Koordinatenprojektionssysteme (CRS) verwendet, um die Erdoberfläche genähert abzubilden, sodass räumliche Entfernungen genauer dargestellt werden können.
Note that our ‘weighted distance’ has just provided a relative difference between the flow accumulation cells and non-flow-accumulation cells. We still must convert these values to some time-dependent form.
For this exercise, we will assume the average velocity of water is 1 m/s value for the flow accumulation cells, and 0.01 m/s for hillslope cells. Let’s calculate the unweighted and weighted time of concentration (time to travel from furthest cell to the outlet).
The DEM is from the USGS 3DEP program and is roughly 30m resolution (each pixel represents a square area of roughly 30m by 30m). We have reduced the (distance) weight of channel cells, so the distance is proportional to our estimate of hillslope velocity, and we can convert the distance to time by dividing the weighted distance by 1 m/s.
Beachtet, dass unsere „gewichteter Distanz“ nur einen relativen Unterschied zwischen den “Abflussakkumulationszellen” und den “Nichtflussakkumulationszellen” liefert. Wir müssen diese Werte noch in eine zeitabhängige Form umwandeln.
Vorerst gehen wir davon aus, dass die durchschnittliche Fließgeschwindigkeit für die Abflussakkumulationszellen 1 m/s und für die Hangzellen 0,01 m/s beträgt. So können wir die ungewichtete und gewichtete Konzentrationszeit (Zeit von der entferntesten Zelle bis zum Auslass).
Das DEM stammt aus dem USGS 3DEP-Programm und hat eine Auflösung von etwa 30 m (jedes Pixel stellt eine quadratische Fläche von etwa 30 x 30 m dar). Wir haben das (Entfernungs-)Gewicht der Gerinnezellen reduziert, sodass die Entfernung proportional zu unserer Schätzung der Hanggeschwindigkeit ist und wir die Entfernung in Zeit umwandeln können, indem wir die gewichtete Entfernung durch 1 m/s dividieren.
# festlegen der Rasterpixelauflösung (m x m)
resolution = (30, 30)
# Die Zellen können nach ihrem gewichteten Abstand zum Auslass gruppiert werden, um den
# den Prozess der Berechnung des Beitrags jeder Zelle zum finalen Abfluss am Auslass zu vereinfachen
dist_df = pd.DataFrame()
dist_df['weighted_dist'] = weighted_dist.flatten()
print(f'Dataframe vor der Aufbereitung {dist_df}')
# entferne Zellen, welche nicht zum Einzugsgebiet gehören
# und runde die Abflusszeit auf die nächste Entfernugs Einheit (Anzahl an Zellen)
dist_df = dist_df[np.isfinite(dist_df['weighted_dist'])].round(0)
print(f'Dataframe nach der Aufbereitung: {dist_df}')
Dataframe vor der Aufbereitung weighted_dist
0 inf
1 inf
2 inf
3 inf
4 inf
... ...
68938 inf
68939 inf
68940 inf
68941 inf
68942 inf
[68943 rows x 1 columns]
Dataframe nach der Aufbereitung: weighted_dist
76 46.0
77 47.0
418 45.0
419 46.0
420 43.0
... ...
68898 64.0
68899 64.0
68900 64.0
68901 64.0
68902 64.0
[36526 rows x 1 columns]
# Ermitteln der Anzahl der Zellen für jeden Abstand
grouped_dists = pd.DataFrame(dist_df.groupby('weighted_dist').size())
grouped_dists.columns = ['num_cells']
# Verteilungen der gewichteten Abstände darstellen
W = np.bincount(weighted_dist[np.isfinite(weighted_dist)].astype(int))
fig, ax = plt.subplots(figsize=(10, 5))
plt.fill_between(np.arange(len(W)), W, 0, edgecolor='seagreen', linewidth=1, facecolor='lightgreen', alpha=0.8)
plt.ylabel(r'Anzahl der Zellen bei Entfernung $x$ vom Auslass', size=14)
plt.xlabel(r'Entfernung vom Auslass (x) [Zellen]', size=14)
plt.title('Weite Funktion W(x)', size=16)
Text(0.5, 1.0, 'Weite Funktion W(x)')
flow_velocity = 0.1 # m/s
w_time = weighted_dist * resolution[0] / flow_velocity / 3600
W_time = np.bincount(w_time[np.isfinite(w_time)].astype(int))
fig, ax = plt.subplots(figsize=(10, 5))
plt.fill_between(np.arange(len(W_time)), W_time, 0, edgecolor='seagreen', linewidth=1, facecolor='lightgreen', alpha=0.8)
plt.ylabel(r'Anzahl der Zellen $x$ Stunden vom Auslass', size=14)
plt.xlabel(r'Zeit zum Auslass $x$ (Stunden)', size=14)
plt.title('Weite Funktion W(x) ~ Einheitsganglinie', size=16)
Text(0.5, 1.0, 'Weite Funktion W(x) ~ Einheitsganglinie')
Hinweis: Oben haben wir eine geschätzte Einheitsganglinie für 1 Niederschlagseinheit entwickelt, die gleichmäßig über das Einzugsgebiet fällt. Dies haben wir allein mit Höhendaten (DEM) und einer Annahme zum Unterschied in den Fließgeschwindigkeiten zwischen Gerinnen und Hangströmung gemacht.
Berechnung des Gesamtabflusses und Vergleich ihn mit dem gemessenen Tagesvolumen für die zweitägige Aufzeichnung#
HINWEIS: Wenn Ihr den Abflusskoeffizienten unten aktualisiert, müsst ihr den Code von hier aus erneut ausführen, um den DataFrame „runoff_df“ neu zu initialisieren. Andernfalls werden die Werte akkumulieren.
# Erstellung einer Einheitsganglinie für jeden Zeitschritt
runoff_df = pd.DataFrame(np.zeros(len(hourly_df)))
runoff_df.columns = ['Total Precip (mm)']
runoff_df.index = hourly_df.index.copy()
end_date = pd.to_datetime(runoff_df.index.values[-1]) + pd.DateOffset(hours=1)
max_distance = max(grouped_dists.index)
def calculate_flow_time(distance, v):
# die Länge des Fließwegs (in Zellen) umrechnen
# in eine Fließzeit
return np.ceil(distance * resolution[0] / v / 3600)
# Konzentrationszeit
max_flow_time = calculate_flow_time(max_distance, flow_velocity)
print(f'Der maximale Abflussweg ist {max_distance} Zellen, was einer maximalen Fließzeit von {max_flow_time} Stunden entspricht.')
extended_df = pd.DataFrame()
extended_df['Total Precip (mm)'] = [0 for e in range(int(max_flow_time) + 1)]
extended_df.index = pd.date_range(end_date, periods=max_flow_time + 1, freq='1H')
# Anhängen der Extra Zeit and den Abfluss dataframe
runoff_df = pd.concat([runoff_df, extended_df])
runoff_df['Runoff (cms)'] = 0
Der maximale Abflussweg ist 102.0 Zellen, was einer maximalen Fließzeit von 9.0 Stunden entspricht.
Die Zelle unten ist eine umfangreiche Berechnung und kann einige Zeit dauern. Was unten passiert, ist, dass wir für jede Niederschlagsstunde den Zeitversatz berechnen, den der Fluss benötigt, um von jeder Zelle zum Auslass zu gelangen. Da nicht jede Zelle die gleiche Zeit benötigt, um zum Auslass zu gelangen, addieren wir den Abfluss jeder Zelle zu einem zukünftigen Zeitpunkt, der der Fließzeit der Zelle entspricht. Wir machen es etwas effizienter, indem wir Zellen mit gleichem Abstand gruppieren.
Nachfolgend nehmen wir einen Abflusskoeffizienten von 0,3 an. D. h. 30 % unseres Niederschlags werden Abflusswirksam sein sein. Später werden wir sehen, wie gut Modellergebnisse mit den gemessenen Werten zusammenpassen, und versuchen, diesen Koeffizienten zu validieren.
runoff_coefficient = 0.3
for ind, row in hourly_df.iterrows():
this_hyd = hourly_df[['precipitation']].copy()
hydrograph = pd.DataFrame()
for weight_dist, num_cells in grouped_dists.iterrows():
weighted_time = calculate_flow_time(weight_dist, flow_velocity)
outlet_time = ind + pd.Timedelta(hours=weighted_time)
# Rundung der Fließzeit auf die nächste Stunde
# Angleichung an die stündlichen Abflussdaten
if weighted_time < 1:
outlet_time = ind + pd.Timedelta(hours=1)
precip_vol = num_cells.values[0] * row['precipitation']
runoff_vol = precip_vol * runoff_coefficient / 1000 * resolution[0]**2
runoff_rate = runoff_vol / 3600 # Umrechnung von m³/h in m³/s
runoff_df.loc[outlet_time, 'Runoff (cms)'] += runoff_rate
runoff_df['day'] = runoff_df.index.day
# Umrechnung von m³/s in m³/h (cmh)
runoff_df['runoff_vol_cmh'] = runoff_df['Runoff (cms)'] * 3600
cumulative_vol = runoff_df[['runoff_vol_cmh', 'day']].groupby('day').sum()
cumulative_vol
| runoff_vol_cmh | |
|---|---|
| day | |
| 28 | 102690.396 |
| 29 | 293762.808 |
# Abflusszeitreihe importieren
whis_flow_df = pd.read_csv(os.path.join(data_path, '08MG026_daily.csv'),
header=1, parse_dates=['Date'], index_col=['Date'])
# whis_flow_df = whis_flow_df[whis_flow_df['2005-09-27': '2005-09-31']]
# whis_flow_df
whis_flow_df = whis_flow_df['2005-09-27':'2005-09-30'].round(1)
whis_flow_df
| ID | PARAM | Value | SYM | |
|---|---|---|---|---|
| Date | ||||
| 2005-09-27 | 08MG026 | 1 | 1.3 | NaN |
| 2005-09-28 | 08MG026 | 1 | 1.5 | NaN |
| 2005-09-29 | 08MG026 | 1 | 4.8 | NaN |
| 2005-09-30 | 08MG026 | 1 | 2.8 | NaN |
cumulative_vol.values.flatten()
whis_flow_df.loc['2005-09-28':'2005-09-29']
| ID | PARAM | Value | SYM | |
|---|---|---|---|---|
| Date | ||||
| 2005-09-28 | 08MG026 | 1 | 1.5 | NaN |
| 2005-09-29 | 08MG026 | 1 | 4.8 | NaN |
# für den Abfluss in Whistler kann der Basisabfluss subtrahiert werden
# nehmt an, dieser ist 1.3 cms für 27st September
whis_flow_df['excess_flow'] = whis_flow_df['Value'] - 1.3
# Umrechnung des durchschnittlichen täglichen Durchflusses in das tägliche Volumen in m³/Tag
whis_flow_df['excess_volume'] = whis_flow_df['excess_flow'] * 3600 * 24
whis_flow_df.loc['2005-09-28':'2005-09-29', 'model_volume'] = cumulative_vol.values.flatten()
fig, ax = plt.subplots(1, 1, figsize=(16,4))
ax.plot(runoff_df.index, runoff_df['Runoff (cms)'], label="Modellierter stündlicher Abfluss")
# ax.plot(whis_flow_df.index, whis_flow_df['Value'], label='Gemessener täglicher Durchschnittsabfluss [m³/s]')
ax.set_xlabel('Datum')
ax.set_ylabel('Abfluss [m³/s]', color='blue')
ax.set_title('Beispiel Niederschlags-Abfluss-Modell')
ax.legend(loc='upper left')
ax.tick_params(axis='y', colors='blue')
ax1 = ax.twinx()
# darstellen der Niederschlagsdaten
ax1.bar(hourly_df.index, width=pd.Timedelta(hours=1),
height=hourly_df['precipitation'], color='green', alpha=0.5,
label="Täglicher Niederschlag")
ax1.set_ylabel('Niederschlag [mm]', color='green')
ax1.tick_params(axis='y', colors='green')
ax1.legend(loc='upper right')
<matplotlib.legend.Legend at 0x10e0747ffa0>
Vergleich des modellierten vs. gemessenen Abfluss#
Wir können den stündlich modellierten und den täglich gemessenen Durchfluss nicht direkt vergleichen, aber wir können das Gesamtvolumen vergleichen.
Wir vergleichen also den modellierten Überschussabfluss aus Niederschlägen mit dem gemessenen Abfluss (bereinigt um den Basisabfluss). Nachfolgend subtrahieren wir den Basisabfluss (der als erster Tag des gemessenen Abflusses angenommen wird) von der Zeitreihe des täglichen Abflusses, addieren den gesamten Volumenfluss über das zweitägige Niederschlagsereignis und vergleichen ihn mit dem gesamten modellierten Abflussvolumen (der bereits vorhanden ist) des selben Zeitraumes.
tot_runoff = whis_flow_df[['excess_volume', 'model_volume']].sum()
relative_error = 100*( tot_runoff['excess_volume'] - tot_runoff['model_volume']) / tot_runoff['excess_volume']
print(f'Bei einem Abflusskoeffizienten von {runoff_coefficient} ({100*runoff_coefficient}% des Niederschlages sind effektiver Niederschlag, der Vorhersageerror über das zweitägige Niederschlagsereignise ist {relative_error:.0f}%) ')
Bei einem Abflusskoeffizienten von 0.3 (30.0% des Niederschlages sind effektiver Niederschlag, der Vorhersageerror über das zweitägige Niederschlagsereignise ist 12%)
Ähnlich bei bei dem Zeitbeiwertverfahren schauen wir uns die Unsicherheit an, welche durch die Wahl des Abflusskoeffizienten begründet ist. Die SCS-Verfahren drücken Niederschlag, der nicht effektiv ist (also nicht als Direktabfluss wirksam wird), als Verlust aus.
In diesem Beispiel testen wir den Niederschlagsverlust durch Infiltration, der den effektiven Niederschlag, die Verzögerungszeit und die Konzentrationszeit bestimmt. Um die Annahmen besser zu veranschaulichen, visualisieren wir den Niederschlag in einer Bereich von (konstanten) Verlustraten im Vergleich zu den stündlichen Daten, die wir oben importiert haben.
precip_losses = [1, 2, 4] # mm/hr
p = figure(title=f'Niederschlag vs. Verluste', width=750, height=300, x_axis_type='datetime')
p.vbar(x='time', width=pd.Timedelta(hours=1), top='precipitation',
bottom=0, source=hourly_source, legend_label='Stündlicher Niederschlag',
color='royalblue', fill_alpha=0.5)
t_loss = [hourly_df.index[0], hourly_df.index[-1]]
colors = ['blue', 'green', 'orange', 'red']
i = 0
for l in precip_losses:
p.line(
t_loss, [l, l],
line_dash='dashed', legend_label=f'{l}mm/h Verlust',
color=colors[i]
)
i += 1
p.legend.location = 'top_left'
p.xaxis.axis_label = 'Datum'
p.yaxis.axis_label = 'Niederschlag [mm]'
p.toolbar_location='above'
show(p)
Das obige Diagramm zeigt den effektiven Niederschlag als blauen Bereich über jeder der Verlustfunktionslinien. Wir können die verschiedenen Mengen des effektiven Niederschlags jedoch auch grafisch darstellen.
for l in precip_losses:
hourly_df[f'excess_{l}'] = (hourly_df['precipitation'] - l).clip(0, None)
hourly_df['day'] = hourly_df.index.day
excess_rainfall = hourly_df[[f'excess_{l}' for l in precip_losses] + ['day']].groupby('day').sum()
excess_rainfall
| excess_1 | excess_2 | excess_4 | |
|---|---|---|---|
| day | |||
| 28 | 8.0 | 1.8 | 0.0 |
| 29 | 10.4 | 4.6 | 1.4 |
hourly_source = ColumnDataSource(hourly_df)
p = figure(title=f'Effektiver Niederschlag', width=650, height=250, x_axis_type='datetime')
i = 0
for l in precip_losses:
p.vbar(x='time', width=pd.Timedelta(minutes=20), top=f'excess_{l}',
bottom=0, source=hourly_source, legend_label=f'{l}mm Verlust',
color=colors[i], fill_alpha=0.5)
# Verschiebung des Index um 20min, damit die Balken nebeneinander liegen
hourly_df.index += pd.Timedelta(minutes=20)
hourly_source = ColumnDataSource(hourly_df)
i += 1
# in jedem Fall den Index zurücksetzen
hourly_df.index -= pd.Timedelta(hours=1)
p.legend.location = 'top_left'
show(p)
Finally, let’s check the range of precipitation loss values against the total precipitaton volumes, as well as the measured total runoff volume, and see how these compare to our assumed runoff coefficient of 30%.
Lasst uns abschließend den Bereich der Niederschlagsverlustwerte anhand der Gesamtniederschlagsmengen sowie der gemessenen Gesamtabflussmenge vergleichen um zu sehen, wie diese mit unserem angenommenen Abflusskoeffizienten von 30 % übereinstimmen.
total_hourly_vol = hourly_df.groupby('day').sum()[['precipitation', 'excess_1', 'excess_2', 'excess_4']]
tot_df = total_hourly_vol.sum()
for ev in [1, 2, 4]:
e_pct = tot_df[f'excess_{ev}'] / tot_df['precipitation']
tot_df[f'pct_excess_{ev}'] = e_pct
print(f'{ev}mm Abflussverlust führt zu {100*e_pct:.0f}% des Niederschlags, der als Direktabfluss wirksam wird')
1mm Abflussverlust führt zu 46% des Niederschlags, der als Direktabfluss wirksam wird
2mm Abflussverlust führt zu 16% des Niederschlags, der als Direktabfluss wirksam wird
4mm Abflussverlust führt zu 3% des Niederschlags, der als Direktabfluss wirksam wird
Fragen, die im Rahmen der Hausarbeit beantwortet werden sollen:#
Welche Faktoren in unserem Einheitsganglinienmodell wirken sich auf die Größe des Scheitelabflusses aus?
Wie lassen sich die Schätzungen des SCS-Scheitelabflusses mit der Gerinnekapazität vergleichen, die wir zu Beginn ermittelt haben? Würdet ihr das Risiko eingehen, dass das Gerinne überläuft, möglicherweise das aquatische Ökosystem in Fitzsimmons Creek gefährdet und ihr somit das Fischereigesetz verletzen würdet?
Gegen Ende der des Notebooks schätzten wir die Abflussreaktion anhand der täglichen Niederschläge und gehen von einem Abflusskoeffizienten von 0,3 aus. Wie könntet ihr Beweise finden, die diese Zahl oder eine andere Zahl stützen?
Würdet ihr eines der in diesem Notebook besprochenen Modelle als „besser“ oder „schlechter“ hervorheben? Können Sie das beste und das schlechteste Modell in ein paar kurzen Punkten zu folgenden Punkten verteidigen: i) Informationsanforderungen, ii) Komplexität, iii) Interpretierbarkeit?
Bonus#
Wir haben die Python-Bibliothek verwendet, um ein Einzugsgebiet an einem bestimmten Ort am FitzSimmons Creek in Whistler abzugrenzen. Angenommen es ist gegeben, dass die Auflösung etwa 30 m betrug, wie würdet ihr die Entwässerungsfläche a)anhand des DEM und b) anhand des Abflussakkumulationsdiagramms ermitteln?
Mega-Bonus#
Der Anteil des effektiven Niederschlags, hängt von vielen Bedingungen ab. Im Notizbuch 1-Datenordner findet ihr Langzeitaufzeichnungen des täglichen Niederschlags und des täglichen Abflusses. Könnt ihr anhand beider Aufzeichnungen einen “langfristigen” (robusteren) Abflusskoeffizienten schätzen? Könnt ihr außerdem einen Bereich von Abflusskoeffizienten festlegen?